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A method of delecting and contct- 
ing iu»Migid body motion in a sequence 
of miages. for instance MR! images of tlie 
human breast TTje method uses a similar- 
ity measure, such as mutual inibnnatiofL 
to estimate the probabilities of a plunl- 
ity of candidate movements for each of a 
phiraUty of sampling polhta in the hnage. 
The probabilities of the candidate move- 
ments are refined in an Iterative process by 
multiplying diem widi wdghled probabO- 
ities of the most probable motions for die 
neighbouring sampling points. After iler. 
adon the modon field is generated by tak- 
ing the movement of the sampling point 
die candidate movement witfi die highest 
probabUiiy after die iteratioa process. Tlie 
sequence of images can be coireci&d by 
die motion field and dien die process re- 
pealed using different, f« instance more 
closely spaced, sampling points for ftmher 
rcfincmcnL The process is particulariy ad- 
vantageous for detecting and conecdng for 
non-rigid movements in Images which do 
not contain recognisable geometric features 
and in images which are non-«onservative 
ix. the total amount of brightness in die 
image changes widi time, for instance as a 
result of the introduction of contrast agent 
and Its dynamic talc»-iip by die tissue be- 
ing imaged. 
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mmOV AND APPARATTT*; FOP T|Vf ^ GE PTincp^t^ r^ 

This inventioa relates to a method and apparatus for processing image data, 
and in particular for processing a sequence of images of a non-rigid body to coirect 
for movement of the non-rigid body between the images. 

There are many applications in which it is useful to detect coircsponding 
features in successive images of a body, for instance to allow detection of the 
movement of the body and to coirect for movement of it Many methods have been 
proposed for detecting such coiresponding features. The successful detection of 
corresponding feattires m successive images allows the registration of the images, i.e. 
allows subsequent images to be transformed to elimmate differences caused by 
motion. 

A common way to detect motion in successive images is to use landmaric and 
segmentation based algorithms (also known as spaise matching). In such raediods, 
geometric features which can be recognised in each image are matched, and these 
features are used to generate a transforaiation which can be applied to the whole 

image. However, this method is not useful for images which lack recognisable 
geometric features. 

For instance, in the field of diagnosis of breast cancer, an imaging technique 
which is cuirendy being explored is magneUc resonance imaging (MRI). Figure 1 of 
the accompanying drawings shows schematically the arrangement for MRI of a 
human breast in which the patient 1 Ues prone inside the main MR coil 3 on a couch 
5. With this arrangemem the breasts hang pendulously inside the double breast coU 7 
as shown diagrammatically in Figure 2. This technique has certain advantages as 
compared with, for instance, x-ray mammography, in that its use is not restricted to 
post-menopausal women or to breasts which do not include scar tissue, and it is also 
a 3-D technique (acmally producing a set of vertically displaced 2-D slices). 
However, it does have certain problems. Firstly, abnormalities in the breast camiot 
be distinguished in a standard MRI image. To detect abnormalities it is necessary to 
inject a paramagnetic compound (typically CJd-DTPA) as a contrast agent and then 
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scan every minute for approximately 7 mimites to monitor the take up of the contrast 
agent. Each scan consists of 16 slice images iii. for instance, the X, Y plane as 
showninFigure 1. Because ofthe length oftime over which imaging occurs, it is ' 
very likely that the patient will move during that time. These movements make it 
i difficult, or sometimes impossible, to interpret the succession of images because not 
only are the images changing through take up ofthe contrast agent, but also because 
of movement ofthe patient For instance Figure 3(A) illustrates a pre-ccntrast image 
and Figure 3(B) a post-contrast image in which the patient relaxed the pectoral 
muscle between the two images. The pectoral muscle can be seen as the darker grey 
tissue towards die top ofthe image. It can be seen that Ae effea of relaxation is to 
change the shape of die breast and cause non-rigid movement m the tissue 
surrounding the muscle. Such non-rigid movement is difficult to compensate for. 

It will also be appreciated fiom Figure 3 that the images are not of the type in 
which there are clearly recognisable geometric features which can help to detect die 
movement between images. Further, because die contrast agent is being taken up 
during the imaging process, features which have one level of brightness in a first 
image will not necessarily have die same value of brightness in a subsequent one (die 
image is said to be non-conservative). In fact, the features which are of most interest 
(abnormalities such as tumours) take up the contrast agent quicker dian odier tissue 
and so inevitably, and deliberately. wiU be changing in intensity between die 
different images, 

Thus diese images are ofthe type which require correction for motion, but 
which are not susceptible to normal motion detection techniques. 

The present invention provides a technique for detecting and correcting 
successive images for movement where die movement is non-rigid-body movement 
and where tiie images may be tiiemselves changing as a function oftime. It is 
dierefore applicable to images which have non-conservative flow, i.e. die total 
amomit of brightness in die image changes. The invention is not only applicable to " 
MRI techniques, but odier non-conservative image flows where tiiere is movement 
between images. 

The present invention achieves diis by calculating and storing for each of a 
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plurality of sampling points in the image a probabiUty map for the displacement of 
that point between the image and a time separated image. Thus not only one possible 
displacement is considered for each sampling point, but instead all are considered. 
This probability map. namely the probabUities that the sampling point has moved 
5 a particular direction, is then refined in an iterative process by considering the 
probability maps for sampling points neighbouring It 

Thus whereas with previously proposed techniques only one "optimal" 
displacement vector was considered for each sampling pomt. with the present 
invention all are considered and stored for future processing. With the previous 
10 techniques a first estimate of the motion field was made based on the single optimal 
vector for each point, and then the field was refined by iteratively recomputing those 
vectors. However, in non-rigid flow there can be several plausible motions locaUy 
and so algorithms which ignore that fact can fall into local minima which do not 

correspondtotheultimatelyoptimalsolutionduringtheiterativeprocess. Because 
the present invention retains die fact that there are several possible alternatives to the 
"optimal" displacement vector, it can avoid being trapped in incorrect local minima. 

Thus in more detail, the present invention provides a method of processing 
image data of a plurality of time^eparated images of a non-rigid body to detect 
movement of tiie body, comprising die steps of:- 

for each of a pluraUty of sampling points in each image calculating and 
storing a pluraUty of candidate movements together with the estimated probability of 
each candidate; 

iteratively recalculating for each sampling point die probability of each of the 
candidate movement based on the stored probability of that candidate movement and 
the probabilities of the candidate movements at other sampling points; and 

generating from the recalculated probabilities a motion field indicative of tiie 
non-rigid body movement. 

The sampling points are set as a regular array over the image and the said 
other sampling points which are considered in die recalculation process are 
preferably neighbouring sampling points. The sampling points need not correspond 
to individual pixels in the image, but are preferably spaced by several pixels. 
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The estimated probabilities for the candidate movements are found by using a 
similarity measure which indicates the similarity of each sampling point to sampling 
points in the preceding image. Hie simUarity measures can be converted into 
probabilities by simply normalising them for each sampling point so that they sum to 
unity. The candidate movements are, of course, the vector displacements which map 
the respective sampling points in the preceding image to the current sampling point. 

As mentioned above a problem with non-conservative images such as MRI 
images when a contrast agent is being dynamically taken-up, is that points of a 
certain intensity in one image cannot be expected to have the same intensity in 
another. The present invention overcomes this problem by careful- choice of 
similarity measure, and in particular by selecting the similarity measure from mutual 
infonnation, normalised mutual inforaiation or entropy correlation coefficient 

The iterative recalculation of the stored probabilities comprises multiplying 
each stored probability by die product of the stored probabilities for the neighbouring 
sampling points. Preferably the only maximum probability for each of the 
neighbouring pomts is taken, more preferably weighted according to the difference 
between the candidate movement and the maximum probability movement of the 
neighbouring sampling point. The recalculation can also be limited to use only 
movements which are judged to be similar to the movement of the sampling point 
under consideration, for instance where the magnitude of the displacement caused by 
the movements differ by less than a preset amount 

The number of iterations in the recalculation can be set as desired, 
conveniently according to the distance between salient points in the image. 

Once the probabilities of the candidate movements have been refined by the 
iterative process, the motion field is generated by selecting as the movement at each 
sampling point tiiat candidate movement which has die maximum probability. This 
motion field can tiien be used to correct later images in die sequence by converting 
the motion field into a transformation which can be applied to each pixel. The 
motion field can be converted to a transformation eitiier by interpolation or by fitting 
a parametric transformation to the motion field. 

Having corrected the images it is possible to repeat die process using 
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differenUy spaced sampling points, for instance more closely spaced sampling points, 
and preferably to do this iteratively to further refine the detection of the motion field 
and registration of the different images. 

The invention also provides a corresponding apparatus for processing image 
5 data and the invention can be embodied as a computer program which may. for 
instance, be stored on a computer readable storage medium. 

The invention will be further described by non-limitative example with 
reference to the accompanying drawings in which: 

Figure 1 is a diagrammatic representation of the arrangement for MRI 
0 examination of a patient; 

Figure 2 illustrates a typical magnetic resonance image of both human breasts 
and a diagrammatic representation of the position of the double breast coil of the 
apparatus; 

Figures 3(A) and (B) illustrate example pre-contrast and post-contrast images 
of a human breas^ 

Figures 4(A) to (F) illustrate comparisons of diflFerent similarity measures 
under contrast enhancement; 

Figure 5 illustrates a typical MRI image and the intensity profile as aflFected 
by the bias field of the MRI apparatus; 

Figure 6 iUustrates results of the appUcation of an ei^bodiment of the present 
invention to pre and post-contrast MR images; 

Figure 7 illustrates further results of an embodiment of the present invention 
as applied to pre and post-contrast MR images; 

Figure 8 illustrates results of applying an embodiment of the present 
invention to pre and post-contrast MR images; 

Figure 9 illustrate pre and post-contrast images using different similarity 
measures; 

Figure 10 illustrates two images with a mutual displacement; and 

Figure 1 1 shows the similarity scores for two sample points in the Figure 10 

image. 

Any image registration algorithm, i.e. an algorithm for detecting and 
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correaing for movement between images, requires a metiiod of detennining whether 
two sub-images represent the same part of the imaged object In the case of breast 
MRI. this amounts to detennining whether the two sub-images represent the same ' 
anatomical location. As mentioned above, while many prior techniques lely on 
recognising geometrical features, this is not suitable for all images, especially images 
from breast MRI. Magnetic resonance signal intensity is a subject to change between 
images because of the deliberately induced contrast enhancement, and this change 
can be of the order of 100%. Further, die deformation which results from the 
movement is non-rigid. Thus die method of detennining whether the two sub-images 
represent the same physical location is important. In this embodiment of the 
invention this is achieved by using a "similarity measure", of which there are several 
types. 

One measure which can be used is the nonnalised "centred" cross- 
correlation:- 

^ r where z/a — > I(r.) 

This is a measure of die least squares residual error from fitting a line to a set 
of data and consequently for registration can be expect to peifonn best when tiieie 
exists a linear relationship between die image intensities. 

Anotiier similarity measure, which is used in die embodiments of tiie present 
invention because of its suitability for breast MR images, is die "mutual 
infonnation*. Mutual infonnation Correlative entixipy) is astatistical measure tiiat 
describes how likely it is tiiat a random variable 7 is fimctionally related to A! It is 
defmed in tenns of tiie entropy and joint entropy of die random variables Xznd Y. 
Entropy can be interpreted as either the degree of disorder or die mfomtation content 
of jr. It is defined by: 



H{X) = -[p{x)\npix)dx 
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which can be estimated discretely for intensity quanta/?, by: 

i-A'-l 



Just as joint and conditional probability distributions can be used to relate co- 
occurrences of two random variables, joint and conditional entropies relate the 
predictabilities of two random variables. Conditional and jomt entropies ai« defined 
respectively as: 

ffiY,X) = -ljpiy,x)lnp(y,x)dxefy 

A'-l A'-l 



Conditional entropy is a measure of the randomness of Kgiven knowledge of 
a: AsKbecomes more dependenton^/friW gets smaller. On its own. this is not 
a measure of dependency, since a small value of HOVQ may only imply Hfrj is 



small. 



The mutual information is defined as: 

MiX,Y)^H{Y)-H{Y\X) 

The conditional entropy can be expressed in terms of the joint entropy: 
H(Y\X)=H{XJ)-H{X) 



SO that the mutual information can be written as: 
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M{XJ)^H{X)^ H(Yh H{XJ) (I) 

In maximizing this measure, transformations are encouraged such that the two 
sample windows are functionally related. 

Mutual infonmation is consequently a statistical measure for which the only 
assumption required is that two pixels that have the same intensity in the first image 
also have the same intensity in the second image. In Breast MR images, this is not 
always true over the whole image - places where it does not hold are of diagnostic 
interest. It is, however, likely to be true on a local scale over most of the unage. 
Thus this is a measure which, over most of the image is expected to be robust under 
dynamic contrast enhancement making it highly suitable for MRI. 

A variation of this is the measure of nonnalised mutual information M„(X,Y). 
which is defined as the ratio of joint and marginal entropies of the images and is 
designed to consider the amount of mutual infonnation with respect of the 
infonnation provided by the individual images:- 

Without this nonnalization, a registration algorithm will favour transfonnations 
which result in high marginal entropies of the two images, rather than 
transfonnations which maximize tine proportion of infonnation shared between the 
images. The only problem with the nonnalization is that it results in two image 
locations which are similar but have very little infonnation (i.e. unifonn intensity) 
being considered as important as two locations which are rich in information (i.e, 
textured). 

A further variation on mutual infonnation which again seeks to normalise 
with respect to the infonnation in die individual images is the entropy cowelation 
codficient which is defined as:- 
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The present invention can use any of these as the similarity measure, the choice being 
dictated by the nature of the images being corrected. 

Figure 4 illustrate the results of applymg these different types of similarity ' 
measure to a test image. The two test images are shown in Figure 4. Image 1 has an 
intensity gradient in the x-direction and two artificial features (to simulate ducts or 
blood vessels). Image 2 was created from image 1 by translating it by 5 pixels, 
rotating it by 10' and increasing the intensity by 40% (the 2 features were different 
enhancements - 100% and 10% respectively) and the intensity gradient present in the 
fust image was reversed in direction. The similarity measures were tested by 
selecting a point for the similarity kernel in image 2 and translating the kernel in 
image I by an amount -20<dx< 20 and calculating the similarity measures between 
the two kernels. The correct result is a maximum similarity at A = +5 - the 
translation between the two images. 

Figure 4 also shows the results of this experiment. Note that with normalised 
15 cross-correlation, only positive correlations (enhancement) are of interest so the 

graph shows die rectified (minimum vahie is zero) cioss-coirelation. As can be seen 
fiora the variations in the computed similarity measures, all four measures : 
normalised cross-coirelation, mutual information, normalised mutual information 
and die entropy correlation coefficient have a maximum at A = 5 (the correct 
displacement between the images). On the basis of these results alone it is difficult 
to choose between the similarity measures for use in contrast enhanced breast MR 
images (and indeed in tests on real images, good results have been achieved with all 
four similarity measures). Howeverthechoiceof similarity measures is made on the • 
basis ofunderstanding what each ofthe similarity measures mean. The difference 
25 between these measures is that normalised cross-correlation measures the error or 
disagreement between two locations whereas the three information measures 
compute the supporting evidence or agreement. Contrast enhancement causes local 
areas of disagreement (changes in the image) where the contrast is taken up strongly 
i.e. the area is potentially suspicious. A registration algorithm based on normalised 
cross-correlation would attempt to reduce diese areas of disagreement, potentially 
causing a tumour to be missed. Consequently, mutual information is a more suitable 
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similarity measure for breast MR image registration. For otlier applications than . 
breast MRI other similarity measures are usable. For instance, nonnalised cross 
correlation can be used for ultrasound images or standard (ie CCD) images. 

Having decided on the similarity measure to be used, it is tiien necessary to 
use it to detect the non-rigid motion field. 

In the following, the similarity measure (chosen from those mentioned above) 
between two patches at {x^yYin image 1 and (x-^ti^y^vY in image 2 will be expressed 
as S(r,u) where r - (x^/znd u = («, vY. 

Firstly the two images which are to be compared and registered (ie the second 
will be transformed to eliminate the effects of motion) are processed using the 
selected similarity measure to produce a set of possible (or "candidate") movements 
for each of a plurality of sampling points across the image together with the 
probabilities of each of those candidates. The way this is achieved using mutual 
infonnation is that firstly a desired grid of sampling points is selected. For instance, 
in the examples described later on and illustrated in Figures 6 to 9, at a first scale the 
sampling points were selected as every third pixel. Then for each of these sampling 
points in the fint image the mutual information for the possible displacements of that 
point are calculated. This involves the calculation of the tiiree terms H(X), H(Y) and 
HPC,Y) of equation (1) by summing over a certain number of pixels (the "kernel" 
size) around the sampling point the product of the image intensity with its natural 
log. The term HOO (the entropy in the first image) is. obviously, a constant for a 
given sampling point, while the term H(r) (the entropy in the second image) varies 
for each possible displacement being considered. In the examples shown in Figures 
6 to 9 the kernel size for calculating the mutual information was taken as 15 x 15 
pixels. Thus the product of the image intensities with their natural log is summed 
over a 15 X 15 pixel window centered on the sampling point m the fint image, and 
the same is done for respective 15x15 pixel windows centered on respective 
differently displaced pixels in the second image. The number of different 
displacements considered and the scale is chosen according to the expected motion. 
In the examples shown in Figures 6 to 9 displacements of up to 6 steps (one step is 
the distance between sampling points) in any direction were considered (this is 
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enough because the movement of the breast is constrained by the double breast coil). 
Thus, for instance, in the second image the mutual information of respective 15x15 
windows centered on the pixels displaced by steps of (-6. -6), (-5, -6), (-4. -6) etc. ' 
from the position of the sampling point in the fust image are calculated. An example 
5 of this for two sampling points on a simple image is illustrated in Figures 10 and 1 1. 
In Figure 10 image 1 is a square \vhich is. in image 2, translated slightly to the right 
and slightiy vertically. Figure 1 1 illustrates the similarity scores for the possible 
displacements of two sampling points A and B in image 1 . For point A, the top 
, lefthand comer of the square in image I, the similarity score for a displacement of (3, 
5 3) is maximum (i.e. 1 00) and this is unique because point A can be identified 

unambiguously in image 2 as being the top lefthand comer. Forpoint B. on the other 
hand, its position along the bottom edge is ambiguous and there are several possible 
displacements (illustrated by the multi headed anow in Figure 10) each of which has 
a maximum score of 100. These displacements are (0. 3), (1, 3). (2, 3). (3, 3), and (4, 
3). Other displacements have lower scores. 

It will be appreciated that the arrays of scores shown in Figure 1 1 are for just 
the two points A and B in image 1, and in fact it is necessary to calculate the 
corresponding array of scores for each of the sampling points over image 1 . Each of 
the scores in the array corresponds to a possible displacement (coiresponding to its 
position in the array) and its score corresponds to the probability of that displacement 
being the correct one. 

Thus these initial similarity scores are considered as estimates of the "prior" 
probabilities of those motions. Then these estimates are refined in an iterative 
probabilistic way. In this embodiment Bayes' theorem is used locally to update the 
estimated probability of each motion given the current motion field and a model of 
the deformation. Bayes' theorem can be expressed as: 

^(y„).W2£(Q 
P(x) 

where /'(r|x) is the probability that the answer or "state" is K given the 
measured values x. On the right hand side of the equation. /'(x| 7) is the probability 
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of the measurements x given that the "state" is 7 (the class c»nditional probability). 

Bayes' theorem can hmce be used to recompute (or update) the probability of 
a motion ti at z, given the sunrounding motion: 

nnjus^) a /»(u^|uj?(u,) (2) 

where P(u,) is the prior probability of displacement (motion) u at location x 
and us^ denotes the motion flows of the surrounding locations. The denominator 
PCx) of the right hand side of Bayes' Theorem is a normalization factor, which results 
in the sum of the posterior probabilities being I.O. The same thing is achieved in the 
present embodiment by noraializing the estimates ofPfuJ after eac* iteration: 

Thus to update the probability using the equation (2) above a model of 
Pi^fM is required which expresses the model of drformation. Considering only 
models whose conditional dependencies are locally dependent: 

where Ufc denotes the motion flows at locations immediately adjacent to x. 

An assumption is then made that the probabilities of motions at two different 
neighbouring sites given the tnie motion at x are independent of one another. With 
this assumption, the joint conditional probability is the product of the probabilities at 
each neighbouring site: 

n«*i«,)«n^(«,.*iu,) 

The probabilities of the different motions at a single location, on the other 
hand, are not independent because they are all part of the same distribution. Thus in 
the present embodiment taking the maximum one as an estimate of /'(u,^^| u,) leads 
to a suitable model for-P(u5te|iix): 
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wheKFfv^J is the probability of motion v at x+6x and N, is the set of 
motions similar to u (the set of motions such that |u-«|<£'). Thus the model in 
equation (4) considers each neighbouring location x+6x and for each neighbouring 
location that one of the of similar motions u+6u which has the maximum probability 
is taken, weighted by a Gaussian (the exponential term). The probability of the 
surrounding motions given u as the coirect motion is estimated as the product of 

these weighted maxima from each location. 

It will be noted that the equation (4) requires multiplications and taking. 

maxima. It is therefore computationally cheaper to take logarithms of the similarity 

data so that all multiplication opecatora become additions. The model equation (4) 

then becomes: 




and the update equation is given by substituting this in a log version of 
equation (2): 



which gives 

(5) 



After applying equation (5) iteratively. the resulting motion field can be 
computed by taking as the estimate of the motion at each sampled point in the image, 
the motion corresponding to the maximum posterior probability: 
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The complete algorithm can be summarized as follows:- 

1. Compute simUarity scores using a K x K pixel kernel (K = 15 in this 
example) at a set of equally spaced sample locations (grid size = 3 pixels in this 
example) across the image for all possible motions (actually up to a maximum 
displacement of 6 steps of 3 pixels in this example) according to equation (2). 

2. Normalize the set of similarities at each location to estimate the prior 
probabilities P(u^ using equation (3). 

3 . Take natural logarithms of the prior probabilities . 

4. Compute the posterior probability of each motion given its 
surrounding estimates using the update equation (5). 

5. Re-normalize the estimated probabilities. Go to Step 4 unless the 
required number of iterations (20 in the examples shown in Figures 6-9) have been 
completed. The number of iterations is dependent on the distance information must 
propagate to inteipolate between features. 

6. Find the motion field by taking, as the estimate of the optic flow at 
each location, the displacement corresponding to the maximum posterior probability 
using equation (6). 

The types of deformations allowed by the algorithm are implicit in the model 
for -flCuaJuJ. The model is designed to minimize — . The behaviour of the 

Sc 

algorithm is influenced by the following parameters: 

The variance a„ in P(u4juJ controls the probability of in — 

Sc 

The grid sampling frequency (the spacing over which motion estimates are 
Si 

taken) scales Jx in — . Hence, for a given samplmg frequency Sx, the variance a 
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can be computed to control — . 

& 



) 



The similarity scores are computed at a set of equally-spaced locations across 
the images. The number of iterations required to compute the posterior probabilities 
depends on the sample frequent^ of the motion flow estimates. 

Thtis using the above technique the probability for each of the several 
candidate movements at each sampling point is refined in an iterative process based 
on the maximum probable movement of the surrounding sampling points. When 
the iterative process has finished the motion field is taken as being that given by 
the most probable candidate movement of each sampling point. 

It -will be appreciated that the sampling points do not correspond to 
individual pixels. Thus in order to correa an image it is necessary to have an 
evaluation of the movement for each pixel of the image. A simple way of 
calculating the movement for each pixel would be to interpolate the motion 
veaors at the nearest sampling points to the pixel. However, a better medaod is to 
fit a parametric transformation to the motion flow field which gives a smooth, 
continuous transformation which can be evaluated at any location. This method is 
described and developed by Dcclerck in Declerck, J., Subsol, G.. Thirion, J., and 
Ayache, N. Automatic retrieval of anatomical structures in 3D medical images. 
Technical report, INRIA. 1995 and Declerck, J., Subsol, G., Thirion. J.P., and 
Ayache, N. Automatic retrieval of anatomical structures in 3d medical images.. In 
Ayache, N., editor. First international conference on computer vision, virtual reality 
and robotics in medicine. CVRh4ed'95. Nice. France 1995. Sprmger-Verlag. Lecture 
Notes in Computer Science, which are hereby incorporated by reference. 

The registration algorithm can be further improved by repeating the process 
using different, more closely spaced sampling points, hi the examples of Figures 6- 
9, one repeat was made using a grid size (sampling point spacing) of 3 pfacels, a 
maximum displacement of 3 in steps of 2 pixels, the same 15 x 15 kernel and only 10 
iterations. Thus the initial processing at a fust scale above computes an approximate 
transformation by sub-sampling die images and sub-sampling the possible 
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displacements. Subsequent repeat processing uses progressively finer sampling of 
the images and displacements. 

Thus the complete algorithm can be summarised as foUows:- 

1 . Compute similarity scores at a set of locations across the image for 
sampled possible displacements according to equation (2). 

2. Normilize the set of similarities at each location to estimate the prior 
probabilities P{u^ using equation (3). 

3. Take natural logarithms of the prior probabilities. 

4. Compute the posterior probability of each optic flow given its 
surrounding estimates using the update equation (5). 

5. Re-normalize the estimated probabilities. Go to Step 4 unless the 
required number of iterations have been completed. 

6. Find die motion field by taking, as the estiinate of the optic flow at 
each location, the displacement corresponding to the maximum posterior probability 
using equation (6) and compute a parametric transformation between the images. 

7. Transform the subsequent image to obtain the motion corrected 
image /j*^ 

8. Go to step 1 and repeat at the next finer scale. 

In certain cases, MR images are corrupted by substantial "bias" field. Figure 
5(B) illustrates the variation in MR signal level across of region of fat in the breast 
for the scan shown in Figure 5(A). The present embodiment uses a model-based 
technique to correct for the bias field. 

When the breast moves, its position in the bias field B(x,y) changes. 
Consequently, even without enhancement, there is an expected intensity change. In 
the presence of movement, the apparent contrast enhancement is computed as: 

-W(^ + r/,y + v) 



The true enhancement is 
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Iposix^u,y-\^v) B{x,y) 
hr^i^.y) 5(jc + w,j;+v) 

The effective change in Post-Contrast Signal level due to movement through 
the bias field (the signal change in with respect to the difference between the 
apparent and true contrast enhancement) is therefore: 

A/ = Ipos,{x + z/, + v) - Ipcst{x + u,y + v) ^^^^V) . 

B{x + u,y^v) ^ ^ 

= Mjc + z/,;;-hv)|l- ^ ^^^'^^ ^1 /P^ 



For registration of images, the intensity change due to movement through the 
bias field can be considered as an enhancement (although this may be negative). It is 
important to note this effect, however, since the purpose of motion correction is to 
allow accurate computation of the contrast enhancement in the breast and the effect 
of the bias field on the motion-corrected hnages will result m an apparent mcrease or 
decrease in the amount of contrast enhancement Consequentiy, whether or not tije 
images are bias-corrected before computmg the motion field (it makes very littie 
difference to the performance of the algoritimi) diey must be bias-corrected before 
correcting for tiie computed motion field, or at least compensated for the motion 
through the bias field using equation (8). 

Figures 6 to 9 illustrate tiie application of the present invention to actual MRI 
sequences. 



wo 00/57361 



PCt/GBOO/01020 



-18- 

Figure 6 illustrates results from a patient who moved by approximately 10 
mm during the scanning process. Figure 6A illustrates the pre-contrast image and 
Figure 6B the post-contrast image. Figure 6C illustrates the result of subtracting one 
from the other without correction for this motion. In reality the sequence was near 
impossible to interpret because of the movement. Figure 6D illustrates the motion 
field as calculated by the present invention and Figure 6E the corrected post-contrast 
image. The result of subtraction of this corrected image from the pre-contrast image 
is shown in Figure 6F, which is visibly clearer than Figure 6C. The intensity profiles 
shown at the bottom of Figure 6 illustrate the results of the correction graphically. 
These demonstrate, in particular, the accuracy of registration at the left breast edge. 
The horizontal displacement between the resulting breast edges in the corrected 
subtraction image is less than half a pixel. 

In Figure 7, Figures 7A and 7B show the pre-contrast and subtraction images 
respectively. A tumour is present which is circled in the subtraction image. The 
patient moved approximately 2 mm between the pre and post-contrast images and 
Figure 7C illustrates the corrected subtraction image. The image surfaces are 
illustrated at the bottom of Figure 7 and the surface is more strongly peaked in the 
conected image, the peak corresponding to the tumour. 

In Figure 8 the invention is applied to a post-contrast image for a patient who 
relaxed the pectoral muscles between the pre and post-contrast acquisitions. Figure 
8A shows a pre-contrast image, Figure 8B a post-contrast image and Figure 8C a 
subtraction image without correction. The motion field calculated by the present 
invention is shown in Figure 8B and the correspondingly corrected post-contrast 
image in Figure 8E Subtraction of this corrected image fix)m the pre-contrast image 
results in Figure 8F. In this sequence the movement is extremely non-rigid in that 
the pectoral muscle changes shape and moves approximately by 1 5 mm. There is 
also a tumour which moves about 7 or 8 mm within the breast, but the movement is 
confined almost entirely within the breast, the edge moving about only 3 mm 
(demonstrating the significant non-rigid nature of the problem). Using the invention 
there is only a 1 -pixel error between the edge in the pre and corrected post-contrast 
image as shown in Figure 8F and illustrated in the split images of Figures 8G and 
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8H. Figure 8G shows a pre-contrast/post-contrast split image without motion 
correction, whereas Figure 8H show the corresponding split image with coirection. 
It can be seen that the boundaries of the various features are aligned in the corrected' 
image. 

In the discussion of simikrity measures above, several different possibilities 
were mentioned. Figure 9 iUustrates the results of using the different similarity 
measures in the present invention. Figure 9A shows the basic pre-contrast image and 
Figure 9B the post-contrast image. Figures 9C, D and E show respectively the 
motion field, corrected post-contrast image and corrected subtraction image using 
normalised cross-correlation. Figures 9F, G and H show the motion field, coirected 
post-contrast image and corrected subtraction image respectively using mutual 
inforaiation. It can be seen that in Figure 9E using nomialised cross-correlation, the 
algorithm has computed the motion incorrectly in the region around the tumour, 
which is precisely the location in which it is required to be most accurate. Using 
mutual information, however, no such problems occur. 
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1. A method of processing image data of a plurality of time-separated - 
images of a non-rigid body to detect movement of the body, comprising die steps of- 
fer each of a pluraUty of sampling points in each image calculating and 
storing a plurality of candidate movements together with the estimated probability of 
each candidate; 

iteratively recalculating for each sampling point the probability of each of the 
candidate movement based on die stored probability of diat candidate movement and 
the probabilities of the candidate movements at other sampling points; and 

generating from the recalculated probabilities a motion field indicative of the 
non-rigid body movement. 

2. . A method according to claim 1. wherein the otiier sampling points are 
the neighbouring sampling points. 

3. A method according to claim 1 or 2, wherem the sampling points 
correspond to unit areas containing a set of pixels of the image. 

4. A method according to claim 1, 2 or 3, wherein the plurality of 
candidate movements and their probabilities are calculated by calculating a similarity 
measure indicating the simUarity of each sampling point to sampling points in the 
preceding image, and noraializing the similarity measures to sum to unity, the 
nonnalised similarity measures being stored as said probabilities, and the candidate 
movements are die corresponding vector displacements which map the sampling 
point to the respective sampling points in the preceding image. 

5. A method according to claim 4 wherein the similarity measure is 
selected from mutual information, normalised mumal information, entropy 
correlation coefficient and centered cross correlation. 
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6. A method according to any one of the preceding claims, wherein the 
iterative recalculation of the stored probabilities comprises multiplying each stored 
probability by the product of the stored probabilities for the neighbouring sampling- 
points. 



7. A method according to claim 6 wherein the iterative recalculation of 
the stored probabilities comprises multiplying each stored probability by the product 
of the respective maxima of the stored probabilities for the neighbouring sampling 
points. 

8. A method according to claim 6 wherein the iterative recalculation of 
the stored probabilities comprises multiplying the stored probability for each 
candidate movement at each sampling point by the product of the respective maxima 
of the stored probabilities for the neighbouring samplmg points, the maxima being 
weighted according to the difference between said candidate movement and the 
respective stored movements corresponding to the maxima. 

9. A method according to claim 8 wherein the weighting is Gaussian. 

10. A method according to any one of the preceding claims wherein the 
iterative recalculations of the probabilities for the movements at each sampling point 
use only movements at neighbouring sampling points which are judged to be similar 
to the said movements at each sampling point 

11. A method according to claim 1 0 wherein movements are judged to be 
similar if the difference in magnitudes of the displacements caused by the 
movements is less than a preset amount. 



12. A method according to any one of the preceding claims wherein the 
number of iterations is set according to the distance between salient points in the 
image. 
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13. A method according to any one of the preceding claims wherein the 
motion field is generated by selecting as the movement at each sampling point the 
candidate movement having the maximum probability after said iterative 
recalculation. 



14. A method accordmg to claim 13 further comprising the step of 
correcting the plurality of images for the movement by applying thereto a 
transfonnation based on the motion field, and then repeating the method of claim 13 
using differently spaced sampling points to generate a new motion field. 

15. A method according to claim 14 wherein the transfonnation is 
calculated by fitting a parametric transfomiation to the motion field. 

16. A method according to claim 14 or 15 wherein the steps of correcting 
the plurality of images for the movement and then repeating the method of claim 13 
using differently spaced sampling points, are carried out iteratively with successively 
more closely spaced sampling points. 

1 7. A method according to any one of the preceding claims wherein the 
images are a sequence of magnetic resonance images of a body taking up a contrast 
agent. 

18. A method according to any one of the preceding claims wherein the 
images are a sequence of magnetic resonance images of a human breast taking up a 
contrast agent. 

1 9. Apparatus for processing image data of a plurality of time-separated 
images of a non-rigid body to detect movement of the body, comprising: 

calculating means for calculating for each of a plurality of sampling points in 
each image a plurality of candidate movements together with the estimated 
probability of each candidate; 
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storage means for storing said candidate movements and estimated 
probabilities; 

recalculating means for iteratively recalculating for each sampling point the 
probability of each of the candidate movement based on the stored probability of that 
candidate movement and the probabilities of the candidate movements at other 
sampling points; and 

motion field generating means for generating from the recalculated 
probabilities a motion field indicative of the non-rigid body movement 



20. A computer program storage medium readable by a computer system 
and encoding a computer program for controlling a computer to process image data 
of a plurality of time-separated images of a non-rigid body to detect movement of the 
body by a method comprising the steps of- 
fer each of a plurality of sampling points in each unage calculating and 
storing a plurality of candidate movements together with the estimated probabUity of 
each candidate; 

iteratively recalculating for each sampling point the probabUity of each of the 
candidate movement based on the stored probability of that candidate movement and 
the probabilities of the candidate movements at other sampling points; and 
20 generating ftom the recalcuhited probabUities a motion field indicative of the 

non-rigid body movement. 

21. . A method of processing image data substantially as hereinbefore 
described with reference to and as Ulustrated in the accompanying drawings. 



22. Apparatus for processing image data constracted and arranged to 
operate substantially as hereinbefore described with reference to and as illustrated in 
the accompanying drawings. 

23. A computer program storage medium readable by a computer system 
and encoding a computer program for controlling a computer to process unage data 
by a method substantially as hereinbefore described with reference to and as 
illustrated in the accompanying drawings. 
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(A) Fig.3. (B) 

Pre-contrast Image Post-contrast image 
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(A) 




Fig.8. 
(B) 



Pre-Contrast Image Post-Contrast Image 




(C) 

Subtraction Image 



(D) 



Motion Field 



(E) 



Corrected 
Post-Contrast Image 



(F) 



Corrected 

Subtraction Image 



(G) 

Pre/Post Split Image 




(H) 

Corrected Pre/Post Split Image 
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Fig.9. (A) (B) 

■■'fP!^ifMnii!i[^^^® Post-Contrast Image 
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